#!/usr/bin/perl -w
use strict;
use FindBin;
use Getopt::Long;
use lib "/net/cpp-group/Leo/bin";
use parse_fasta;
use Cwd;
use run_cmd;
use die_error;

#888888888888888888888888888888888888888888888888888888888888888888888888888888

#   Usage

#888888888888888888888888888888888888888888888888888888888888888888888888888888
my $usage = <<"USAGE";

USAGE:

  slr_parallel.pl 

      DESCRIPTION:
              Runs site specific slr analyses in parallel on the cluster. 
  
      REQUIRED PARAMETERS:
              --cdna CDNA.FILE       N.B. These should be the full paths to the files
              --prot PEPTIDE.FILE          
  
      OPTIONAL PARAMETERS:
              --f61                  Use full codon table instead of F3x4
                                          when estimating initial ds tree
                                          N.B. This is more biologically accurate 
                                          but may over-parameterise for short sequences
              --job_name  JOBNAME    Specify this to be able to 
                                          view job progress / kill jobs later


              --msf                  Use alternative formats for peptide
              --clustal                   multiple alignment 
                                          (defaults to FASTA)
              --regex "(.*)"         Regular expression to find 
                                          corresponding peptide and cDNA
                                          identifiers. (See following.)
              --verbose              Print more information
              --working_dir          Perform all analysis here instead of current directory
              --help                 This screen
  

NOTES:
    Run this programme on the cluster head node (currently fgu205 / fgu200).
        (i.e. run "ssh fgu205" to logon etc.)
  

    The corresponding cDNA and peptide identifiers must match.
        By default, identifiers are the first letters until a space/tab.
        Otherwise, the identifiers can be selected using a regular 
        expression (in --regex)

    Always check for errors in the pipeline (All files ending in '.err')
       ***ESPECIALLY***  mapping.err
	   

RESULTS:
    See results.parameters

FAILURES:
    If it fails, it may be because the tree is too large/small.
        Check the pairwise dS values in 'slr.pairwise/pairwise.dS'.
        The corresponding fitch-predicted dS tree is in 'fitch.tree'.

   Happy Hunting.
                  .... Leo (Goodstadt) 1/11/2004 to 17/08/2007

USAGE

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#   Get options

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
# mandatory


# optional parameters
my $help = undef;
my $regex = '(\\S+)';

my $job_name = '';
my $working_dir;
my $prot_file;
my $cdna_file;
my $msf;
my $clustal;
my $f61;
my $verbose;
my $vverbose;
my $queue_name;
my $priority;
{
    open CMDLINE, ">cmd.line" or die "Error:\n\tCould not open cmd.line\n";
    print CMDLINE join (" ", @ARGV), "\n";
    close CMDLINE;
}
GetOptions(
			'p|prot=s'	=> \$prot_file,
			'c|cdna=s'	=> \$cdna_file,
			'regex=s'	=> \$regex,
			'msf'		=> \$msf,
			'clustal'	=> \$clustal,
			'f61'		=> \$f61,
			'verbose'	=> \$verbose,
			'vverbose'	=> \$vverbose,
			'help'		=> \$help,
			'job_name=s'=> \$job_name,
			'working_dir=s'=> \$working_dir,
			'q|queue_name=s'=> \$queue_name,
			'priority=s'=> \$priority,
			);

die $usage if ($help);
die $usage unless ($prot_file && $cdna_file);

my $bin_dir = $FindBin::Bin;
$working_dir ||= getcwd();

die "Error\n\t:[$working_dir/$prot_file] does not exist.\n" unless (-e $prot_file);
die "Error\n\t:[$working_dir/$cdna_file] does not exist.\n" unless (-e $cdna_file);
my $verbose_cmd = $vverbose ? "--verbose" :"";
unless ($regex =~ /^".*"$/)
{
	$regex = '"'.$regex.'"';
}

$f61 = defined ($f61) ? 0 : 3;

$verbose = 1 if ($vverbose);














#_________________________________________________________________________________________
# 
# 	Safe make directory
# 
#_________________________________________________________________________________________
sub my_mkdir($)
{
	my ($dir_name) = @_;
	if (!-e $dir_name)
	{
		mkdir $dir_name
				or die "Error\n\tCould not create the directory [$dir_name]\n$!\n";
	}
	elsif (! -d  "$dir_name")
	{
		die "Error\n\t[$dir_name] is a file not a directory.\n$!\n";
	}

}


sub print_pairwise_ctl;

print "\n\n" if $verbose;


my_mkdir ($working_dir);


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


#		Main logic


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


#_________________________________________________________________________________________
# 
# 	Prepare sequences
# 
#_________________________________________________________________________________________
#
#	Prepare FASTA format protein sequence file
#

sub prepare_fasta_files
{
	if ($msf)
	{
		print
	#STDERR
	 "MSF --> FASTA ...\n" if $verbose;
	run_cmd(<<"CMD");
		$bin_dir/msf_to_fasta.pl $prot_file
			> $working_dir/$prot_file.fa
CMD
	}
	elsif ($clustal)
	{
		print
	#STDERR
	 "CLUSTAL --> FASTA ...\n" if $verbose;
	run_cmd(<<"CMD");
		$bin_dir/clustal_to_fasta.pl $prot_file
			> $working_dir/$prot_file.fa
CMD
	}
	else
	{
	run_cmd(<<"CMD");
		cp $prot_file $working_dir/$prot_file.fa
CMD
	}
	run_cmd(<<"CMD");
		cp $cdna_file $working_dir/$cdna_file.fa
CMD
}




#_________________________________________________________________________________________
# 
# 	Map cDNA onto aa
# 		
# 		Can run locally
# 
#_________________________________________________________________________________________
sub map_prot_to_cdna
{
	print "Using protein alignment to align corresponding cDNA ...\n" if $verbose;
	run_cmd(<<"CMD");
		$bin_dir/map_cdna_onto_aa
			--peptide_sequences $working_dir/$prot_file.fa
			--cdna_sequences $working_dir/$cdna_file.fa
			--peptide_output $working_dir/$prot_file.mapped.fa
			--cdna_output $working_dir/$cdna_file.mapped.fa
			--regex $regex
			--missing_identifiers $working_dir/missing.cdna_identifiers.err
			--err_log $working_dir/mapping.err
			--x_for_aa_mismatches X
			$verbose_cmd
CMD
		unlink ("$prot_file.fa");
	
	if (-f "$working_dir/mapping.err")
	{
		open_or_die(*ERROR, "$working_dir/mapping.err");
		print $_ for (<ERROR>);
	}
	if (-f "$working_dir/missing.cdna_identifiers.err")
	{
		open_or_die(*ERROR, "$working_dir/missing.cdna_identifiers.err");
		print $_ for (<ERROR>);
	}
	
	#
	#	Prepare cDNA onto aa
	#
	{
		print
	#STDERR
	 "Convert mapped cDNA to PAML format ...\n" if $verbose;
	run_cmd(<<"CMD");
		$bin_dir/seq_extract
			--input_file $working_dir/$cdna_file.mapped.fa
			--wrap_line    0
			--output_format p
			$verbose_cmd
			--e $working_dir/extract.err
			> $working_dir/seq.tmp
CMD
	}
	
	if (-f "$working_dir/extract.err")
	{
		open_or_die(*ERROR, "$working_dir/extract.err");
		print $_ for (<ERROR>);
	}
	
}


	#_________________________________________________________________________________________
	# 
	#
	#  -------Create pairwise Ks directory
	#
	#
	#  -------make directories
	#
	#my_mkdir("$working_dir/slr.pairwise");
	my_mkdir("$working_dir/slr");
	
	
	
	
	
	
	#
	#  -------Create PAML control files ..
	#
	print "Create PAML control files ...\n" if $verbose;
	# pre-declaration
#	print_site_specific_ctl();
#	print_pairwise_ctl();
	
	
#sub calculate_pairwise_ds
#{
#	
#	
#	#_________________________________________________________________________________________
#	
#	#_________________________________________________________________________________________
#	#
#	#  -------Calculate pairwise Ks to build tree
#	#
#	
#	
#	print "Calculate pairwise Ks to build tree ...\n" if $verbose;
#	run_queue_cmd(<<"CMD", $job_name."_seq_pairs", $queue_name, $priority);
#	$bin_dir/seq_pairs_kaks
#		-aA
#		-t $working_dir/slr.pairwise
#		-c $working_dir/$cdna_file.mapped.fa
#		--err_log $working_dir/kaks.pairwise.err
#		-p $working_dir/paml.err
#		$verbose_cmd
#		> $working_dir/kaks.results
#CMD
#	
#	if (-f "$working_dir/kaks.pairwise.err")
#	{
#		open_or_die(*ERROR, "$working_dir/kaks.pairwise.err");
#		print $_ for (<ERROR>);
#	}
#
#}
#
#
#sub build_fitch_tree_from_pairwise_ds
#{
#	print "Filter out Ks for build tree ...\n" if $verbose;
#	run_cmd(<<"CMD");
#		$bin_dir/seq_pairs_filter_kaks_results
#			-S
#			-R $working_dir/filter.summary
#			-r 0
#			-f 0
#			-s 0
#			-t 0
#			-n 0
#			-l 0
#			-e $working_dir/filter.err
#			$verbose_cmd
#			< $working_dir/kaks.results
#			> $working_dir/slr.pairwise/pairwise.dS
#CMD
#
#	if (-f "$working_dir/filter.err")
#	{
#		open_or_die(*ERROR, "$working_dir/filter.err");
#		print $_ for (<ERROR>);
#	}
#	
#
#
#	#
#	#  -------Build fitch tree from ks values
#	#
#	print "Build fitch tree from ks values ...\n" if $verbose;
#
#	# run fitch
#	run_queue_cmd(<<"CMD", $job_name."_tree_fitch", $queue_name, $priority);
#	$bin_dir/tree_fitch -i p
#			$verbose_cmd
#			< $working_dir/slr.pairwise/pairwise.dS
#			> $working_dir/fitch.tree
#CMD
#}
#
sub build_tree_with_treebest
{
	# run fitch
#	run_queue_cmd(<<"CMD", $job_name."_tree_fitch", $queue_name, $priority);
	my $tree_best_cmd = <<"CMD";
	/net/cpp-group/cluster/bin/treebest phyml  $working_dir/$cdna_file.mapped.fa > $working_dir/treebest.tree
CMD
	print ($tree_best_cmd . "\n");
	run_cmd($tree_best_cmd);
}	



########################   SLR   ########################

#
#  -------site specific analyses
#
sub site_specific_analyses
{
	print "Site specific analyses with slr ...\n" if $verbose;
	my_mkdir("$working_dir/slr");
	chdir("$working_dir/slr");
	system("$bin_dir/tree_to_slr_format.pl < $working_dir/treebest.tree --sequence $working_dir/seq.tmp > $working_dir/slr/slr.tree");
	system("cp $working_dir/seq.tmp $working_dir/slr/seq.tmp");
	print STDERR "run Slr:\n";

	# /net/cpp-group/tools/slr/bin/Slr
	#print STDERR "$bin_dir/tools/slr/bin/Slr $working_dir/slr/slr.ctl > $working_dir/slr/output.slr\n";
	#run_queue_cmd(<<"CMD", $job_name."_slr", $queue_name, $priority);
	run_cmd(<<"CMD");
	$bin_dir/tools/slr/bin/Slr 
			$working_dir/slr/slr.ctl 
			> $working_dir/slr/output.slr
CMD
	system("mv $working_dir/slr/slr.results $working_dir/slr.results");
	system("mv $working_dir/slr/output.slr $working_dir/slr.output");
	chdir($working_dir);


	print "Analyse results ...\n" if $verbose;

	# summarise SLR output
	{
		open (SLR_STDOUT, "$working_dir/slr.output") or die;
		open (SUMMARY, ">$working_dir/output.summary") or die;
		while (<SLR_STDOUT>)
		{
			last if /^\*\*\*/;
		}
		
		$_ = <SLR_STDOUT>; # ignore next

		# tree
		{
			$_ = <SLR_STDOUT>;
			chomp;
			print SUMMARY "$_\n";
		}
		
		# lnL
		{
			my $lnL = <SLR_STDOUT>;
			$lnL =~ /# lnL = ([0-9\.\+\-\e]+).*/ or die;
			print SUMMARY "$1\tlnL\n";
		}

		# kappa / omega
		{
			my $line = <SLR_STDOUT>;
			$line =~ /Kappa = ([0-9\.\+\-\e]+) Omega = ([0-9\.\+\-\e]+)/ or die;
			print SUMMARY "$1\tkappa\n";
			print SUMMARY "$2\tomega\n";
		}

		# tree parameters
		{
			my $tree = <SLR_STDOUT>;
			$tree =~ /Tree length = ([0-9\.\+\-\e]+), average branch length = ([0-9\.\+\-\e]+) \(min=([0-9\.\+\-\e]+), max=([0-9\.\+\-\e]+)\)/ or die;
			print SUMMARY
			"$1\ttree length\n",
			"$2\taverage branch length\n",
			"$3\tminimum branch length\n",
			"$4\tmaximum branch length\n";
		}
	}

#	run_cmd(<<"CMD");
#	$bin_dir/parse_slr_results.pl
#			--significant_residues  $working_dir/results.significant.fa
#			$verbose_cmd
#			$working_dir/slr.results
#			>> $working_dir/output.summary
#CMD

	

	

	print "All finished...\n\n";

}







#
#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#	parameter files
#







sub print_pairwise_ctl
{
	open PAIR_CTL, ">$working_dir/slr.pairwise/pairwise_kaks.ctl"
			or die "Error\n\t:Could not write to pairwise_kaks.ctl\n$!\n";
	print PAIR_CTL <<"CMD";
seqfile = ../seq.tmp
outfile = results
icode = 0	  			# universal code
noisy = 0
verbose = 0
runmode = -2  			# pairwise
seqtype = 1	  			# codons
CodonFreq = $f61		# F3X4; Could be 3 for codon table Was error fixed thanks to Andreas
aaDist = 0	  			# ??? equal Could be geometric or linear or G1974 Miyata
model = 0	  			# one dN/dS over tree
NSsites = 0	  			# one dN/dS over sequence
Mgene = 0		    	# ??? rates; 1: separate
fix_omega = 0
omega = .4
fix_alpha = 1		 	# gamma distribution shape set to 1 for single rate
alpha = 0		     	# gamma distribution shape set to 1 for single rate
fix_rho = 1
rho = 0
Malpha = 0		     	# ??? different alphas for genes
ncatG = 1		     	# # of categories in dG of NSsites models5D
clock = 0		     	# ??? no clock 1:clock
getSE = 1
RateAncestor = 1	 	# ??? rates alpha > 0 or ancestral rates 1 or 2
Small_Diff = .5e-6
method = 0		     	# simultaneous, 1: one branch at a time
fix_blength = 0		 	# ignore branch lengths from supplied ree
fix_kappa = 0
kappa = 2
cleandata = 0
CMD

}

sub print_site_specific_ctl()
{
	my $file = ">$working_dir/slr/slr.ctl";
	open SLR_CTL, $file
			or die "Error\n\t:Could not write to $file\n$!\n";
print SLR_CTL <<"CMD";
seqfile: seq.tmp
treefile: slr.tree
outfile: slr.results
positive_only: 0
reoptimize: 1
#codonf: $f61
CMD
}


sub print_fitch_parameters
{
	my $file = ">$working_dir/fitch.parameters";
	open FITCH_CTL, $file
			or die "Error\n\t:Could not write to $file\n$!\n";
	print FITCH_CTL <<"CMD";
build_tree.dS
0
7
L
J
33
50
2
Y
CMD
	close FITCH_CTL;
}










#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


#		main


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


prepare_fasta_files();
map_prot_to_cdna();
#calculate_pairwise_ds();
#print_fitch_parameters();
#build_fitch_tree_from_pairwise_ds();
build_tree_with_treebest();
print_site_specific_ctl();
site_specific_analyses();


print "All finished...\n\n";





